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Abstract 

Cells generally change their internal state to adapt to an environmental change, and accord¬ 
ingly evolve in response to the new conditions. This process involves phenotypic changes 
that occur over several different time scales, ranging from faster environmental adaptation 
without a corresponding change in the genomic sequence to slower evolutionary dynamics 
involving genetic mutations and subsequent selection. In this regard, a question arises as to 
whether there are any relationships between such phenotypic changes over the different time 
scales at which adaptive evolution occurs. In this study, we analyzed simulated adaptive evo¬ 
lution in a simple cell model, and found that proportionality between concentration changes 
in adaptation and evolution over all components, and the proportion coefficients were closely 
linked to the change in the growth rate of a cell. Furthermore, we demonstrated that the 
phenotypic variances in component concentrations due to (non-genetic) noise and genomic 
alternations are proportional across all components. These global relationships in cellular 
states were also supported by phenomenological theory and transcriptome analysis of labo¬ 
ratory evolution in Escherichia coli. These hndings provide a basis for the development of 
a quantitative theory of plasticity and robustness, and to determine the general restriction 
of phenotypic changes imposed by evolution. 
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I. INTRODUCTION 


When an environmental condition is changed, biological systems change their state to 
adapt and evolve to the environmental change. Despite the recognized importance to char¬ 
acterize the capacity of adaptation and evolution, discussions on evolvability and plasticity 
have thus far remained at a qualitative, rather than quantitative, level. On the other hand, 
the cellular internal state can now be quantitatively determined by measuring the abun¬ 
dances of a variety of components, including proteins and metabolites. Recent advances 
in high-throughput experimental analysis enable quantification of changes within such a 


Q. 


high-dimensional state space 

After an environmental change, cells may hrst respond by changing the abundances of 
cellular components without changing the genome sequence. The typical time scale of such 
environmental adaptation is generally shorter than several generations. On the other hand, 
over the long-term, i.e., over many generations, the internal state is gradually changed by 
evolutionary dynamics, in which the genome sequence is altered by mutations and individ¬ 
uals with higher fitness are generally selected. Indeed, experimental data of both changes 
in phenotype, reflecting changes in gene expression prohles, and changes in the genomic 
sequence throughout the course of evolution are now available; for exampl^ much data are 
available from results obtained from the experimental evolution of E. coli 2l-l4|. 

Therefore, an important question arises: is there a general relationship between short¬ 
term phenotypic changes in adaptation and long-term phenotypic changes in evolution? Of 
course, the phenotypic changes that occur over different time scales are generally caused by 
different mechanisms, and thus the existence of any relationship between them would be non¬ 
trivial. However, it should be noted that the essence of cellular dynamics is reproduction, 
in which the abundance of each cellular component is roughly doubled, and this constraint 
imposed by cellular reproduction imposes a restriction on the time development in high¬ 
dimensional cellular state space. That is, it is possible that a (glivingh) cellular state would 
be restricted to a sub-space of the high-dimensional state space, described by a relatively 
smaller number of variables. Such restriction to low-dimensional dynamics can provide a 
non-trivial link between the phenotypic changes occurring in adaptation and the long-term 
changes occurring over the course of evolution j^. In fact, some studies have suggested 
that there is a common trend over the thousands of gene expression changes in adaptation 
and evolution, in which genes whose expressions exhibit a larger response to environmental 
change tend to also show a larger response in their expression at the evolutionary scale 
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a, l7|]. Furthermore, global cellular behavior is represented by few macroscopic variables 
such as the growth rate and htness of the system, which govern the entire (high-dimensional) 
dynamics in a cell. Therefore, it is important to uncover the possible relationships in the 
high-dimensional cellular dynamics between the expression of the thousands of proteins and 
metabolites and a macroscopic variable such as the growth rate, throughout the process of 
adaptation and evolution. 

In addition to the phenotypic changes that occur after environmental changes during 
adaptation and evolution, the cellular state also generally exhibits fluctuations even under a 


constant environment and without genomic a. 
tic nature of intra-cellular chemical reactions 


ternations, which originate from the stochas- 
8|, l9| . The possible relationship between such 
non-genetic fluctuations and adaptive responses has also garnered much attention recently. 
The proportionality between such fluctuations and the evolutionary rate of htness and phe¬ 


notype has been demonstrated in bacterial experimental evolution anc 


toy cell models, which are supported by phenomenological theory 12l-ll5|. analogous to the 


in simulations of 


proportionality between the huctuation and response in statistical physics that has been 
well established since Einstein jw, 1I|. In the present context, evolution is considered the 
response to genetic change, and thus the proportionality between huctuation and evolution¬ 
ary rate means that components that are more variable by noise are also more variable by 
genetic change. 

Considering the suggested proportionality between the response of cellular states to the 
environmental change and to genetic change, and also between the response and huctuations, 
one may expect the existence of proportionality among two-by-two quantities, namely, huc¬ 
tuations and responses induced by environmental (non-genetic) perturbations (noise) and by 
genetic changes (mutation). This grand relationship, if conhrmed, is of critical importance 
to evolutionary biology, as it would provide a theoretical basis for the quantitative study on 
the plasticity and robustness underlying adaptive evolution, and could also set a general re¬ 
striction as to the extent of phenotypic changes possible through (future) evolution. Indeed, 
the quantities constituting the relationship can now be measured over high-dimensional cel¬ 
lular dynamics, rehected as changes in the expression of thousands of genes. However, thus 
far, experimental conhrmation of this grand relationship remains premature, and is in need 
of further scrutiny. Accordingly, at this stage, it is important to examine such a relationship 
by adopting an integrative approach combining in silico evolution of a cell model consist¬ 
ing of thousands of chemical species, laboratory evolution of bacteria under environmental 
stress, and phenomenological theory for time development in a high-dimensional state space. 
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In the present stndy, we aimed to uncover such statistical laws underlying the fluctuation 
and response of high-dimensional state variables occurring through adaptive evolution, and 
to connect them with changes in growth rate or htness. 


II. RESULTS 


Evolutionary simulations of a simple cell model 

We employed our previously established mutually catalytic reaction network model, as 
this model is capable of capturing the basic characteristic of cells such as the power-law 
abundances, log-normal fluctuations, fluctuation-response relationship of htness, adaptation 


with fold-change detection, and so forth, in spite of its simplicity j9,ll3|,ll^,ll7|. In the model. 


the cellular state is represented by a set of molecule numbers (W, -^ 2 , • • • , Nk), where W 
is the number of molecules of the chemical species i, which ranges from i = 1 to For the 
internal chemical reaction dynamics, we chose a catalytic network among these K chemical 
species, where each reaction from some chemical i to some other chemical j is catalyzed 
by a third chemical Some resources (nutrients) are supplied from the environment by 
transportation through the cell membrane with the aid of some other chemicals that are 
termed ’transporters’. The environmental condition is given by the concentrations of nutrient 
chemicals. Through catalytic reactions, these nutrients are transformed into cell-component 
chemicals, and a cell divides when the amount of component chemicals reaches a certain 
threshold. Here, to achieve a higher growth rate, the synthesis of the cell components, 
transporters, and chemicals that catalyze the synthesis of those components need to be 
harmonized with the nutrient uptake. We allowed the above toy-cell model consisting of 
catalytic reaction networks to evolve by rewiring the network paths with a given mutation 
rate and selecting the pathways with a certain fraction of cells that showed a higher growth 
rate (See Methods for details). For a given environmental condition, evolution progresses 
so that the cell growth rate, i.e., the inverse of the average division time, is increased 
(Fig. 1). To study the response to environmental change, we then switched the nutrient 
condition after evolution under a hxed condition for 3000 generations (denoted by the arrow 
in Fig. 1). The growth rate initially decreased following this environmental change, and 
then recovered through genetic evolution over generations. Next, we explored the phenotypic 
state changes in response to the environmental and evolutionary changes in order to evaluate 
the relationship between non-genetic and genetic responses upon environmental change. 

As phenotypic state variables for the cell, we computed the abundances of each chemical 
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Ni at the division event. Here, it is convenient to choose X* = log X* as a phenotypic variable, 
since the abundance generally increases exponentially over time through cellular growth, and 
perturbation in a network is also generally amplihed exponentially. Indeed, this choice of 
logarithmic abundances is also relevant to the theoretical argument presented below, as 
well as to transcriptome analysis of gene expression. Note also that the abundances Ni are 
distributed by cells, even for those sharing the same reaction network, due to stochasticity 
in reaction dynamics; thus, the average abundance over all cells is required to study the 
mean response of cells, denoted by (••■). 

After the change in nutrient condition, the abundances of all the components change. 
Let us denote the average change of these abundances by: = (Xj(l)) — (Xj(0)) = 

log i where generation 1 refers to the time point immediately following the envi¬ 
ronmental change, and generation 0 denotes the generation right before this nutrient 
change. Similarly, we dehne the response by genetic evolution after m generations by 
= {Xi{m)) — (Xj(0)). Fig. 2 shows the plot of 6Xf^'" versus for 

m = 5, 10, and 50. Interestingly, proportionality was found between the environmental and 
genetic responses over all components. 

Let us now dehne this proportion coefficient r(m) for across components i. This 

proportion coefficient r(m) is initially close to 1, but with the increase in generations m, it 
decreases towards zero, in conjunction with the recovery of the growth rate. In other words, 
evolution shows a common tendency to abolish the changes in components introduced by the 
environmental change. This common proportionality across all chemicals suggests that the 
proportion coefficient r(m) is a “global variable” over a huge number of chemical species. A 
reasonable candidate for such a global variable is the cell growth rate fi. Hence, it is natural 
to compare the coefficient r(m) with the growth rate. Toward this end, we again computed 
the change in the growth rate = /i(l) — p(0)(< 0) and = /i(m) —/i(0) at the 

m th generation. The ratio gives an index for the recovery in this growth 

rate from the decrease caused by the environmental change, with 0 and 1 representing full 
and null recovery, respectively. In Fig. 3, the proportion coefficient r(m) is plotted against 
this growth rate recovery The proportionality between the two is clearly 

discernible. 

Recalling the possible relationship between fluctuation and response, as is typical in 
statistical physics, we then evaluated whether there also exists a common relationship among 
the variances of all the components. Here, as previously reported j^, the distribution of each 
Ni follows an approximately log-normal distribution, as conhrmed experimentally in the 
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protein abundances in the present cells. Hence, it is again relevant to adopt Xj = log Ni as 
phenotype variable, so that the distribution of Xj follows a roughly Gaussian distribution 
12| . The phenotypic variance Vip{i) for each component i is dehned as the variance of Xj in 
an isogenic population. On the other hand, the variance due to genetic change Vg{i) is dehned 
as the variance of mean Xj over a heterogenic distribution, where the mean is computed 
across clones of a given genotype (i.e., a network), while the heterogenic distribution is 
related to different genotypes (networks) that exist at a given generation. 

The results of simulations are given in Fig. 4, which shows proportionality between V)p(i) 
and Vg{i) across the components i for the evolved population. As the mutation rate is 
increased over evolution, Vg{i) increases as the genotype distribution is broadened, whereas 
Vip{i) remains at the same level, so that the ratio Vg{i)/Vip{i) is increased while maintaining 
the proportionality. This observed proportionality means that the components that are more 
variable owing to noise in the reaction dynamics are also more variable owing to mutation. 

So far, we have conhrmed the existence of common proportionality between non-genetic 
and genetic variances, as well as between the environmental and evolutionary responses. As 
the proportionality between the fluctuation and response is a natural outcome in statistical 
physics, we compared the response and fluctuations in more detail. However, direct com¬ 
parison of the phenotypic variance V)p(i) with the environmental response did not show a 
clearly discernable proportionality. This is probably due to the discrepancy in the dehni- 
tions of the two quantities: the variance originates from very high-dimensional dynamics 
without any specihc directional change, while in the environmental response, only one spe- 
cihc environmental change considering only a few nutrients is applied. To make a more 
direct comparison, we then sampled the environmental responses against a variety of exter¬ 
nal changes introduced by different nutrient conditions to dehne the average environmental 
response = ((5X®”^)^) with (■ ■ ■) over 10^ environmental conditions (see Methods). 

As shown in Fig. SI, this average environmental response showed clear proportionality with 
Vip{i) and Vg{i), respectively. Hence, the proportionality relationships among two-by-two 
quantities, i.e., genetic and non-genetic responses and fluctuations, hold over all components 
(see Fig. 5). 

Thus, the degree of plasticity required to achieve an adaptive response to a new environ¬ 
ment is characterized by fluctuations V)p(i), i.e., those that do not consider environmental 
or genetic changes. On the other hand, when cells are exposed to a novel environment, the 
potential of adaptation is expected to increase. Therefore, when placed in a novel condition, 
it is expected that the phenotypic fluctuations would increase to allow the cells to adapt to 
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the new environment. In Fig. S2, the variances {Vipii),Vg{i)) are plotted before and after 
the environmental change (arrow in Fig. 1). In this case, all of the variances increase while 
roughly maintaining their proportionality. After this increase, the variances decrease over 
generations under a fixed environmental condition, while the proportionality between Vip{i) 
and Vg{i) is maintained. 


Theoretical Argument 

By using a simple cell model, we have conhrmed the common proportionality over thou¬ 
sands of components for genetic and non-genetic variances, and environmental and genetic 
responses. The results suggest the existence of a global variable that governs adaptive evo¬ 
lution. Here, the growth rate /i of a cell is a candidate for such a variable, since, for a cell 
to maintain its composition, every component has to be synthesized in conjunction with the 
growth rate. Indeed, in we considered the dynamics of gene expression 

dxi/dt = fi{{xj}) - jxxi, (1) 


where (iXi gives the dilution of the concentration by the increase in cell volume V, and Xi is 
the concentration of the component i, Xi = Ni/V. By using X, = logXj and Fi{{Xj})xi = 
fi{{xj}), the original stationary state is given by 


F,({x;}) = /i. 


( 2 ) 


Now, with the change in environmental condition E and genetic change G, the expression X* 
is shifted to X* -|- hXj, and /r is shifted to fi + 6fi. Assuming that the change in logarithmic 
concentration is not so large, and taking only the linear part of the changes and using the 
Jacobi matrix Jij = {^)x^=x^, we get 


Jij6Xg{E, G) + yf + yf JG = 6fi{E, G), (3) 

3 

where yf = ^ and yf = respectively. 


Here, G is a coordinate introduced to represent the genetic change. It is not evident 
that the genetic change is represented by only a single variable. However, considering that 
under this scenario evolution progresses under a stressed environmental condition, one could 
project high-dimensional genetic change in the direction required to increase htness (growth 
rate) under the condition, indicating that a single variable G can be introduced; indeed, 
several studies conducted to date support this assumption [^, 13, 151. Accordingly, the 


variable G has the same dimensions as E, and can be scaled so that G and E induce the 






same degree of change in expression. The genetic evolution following the initial stress 5E 
is expected to diminish the environmental stresses, so that evolution occurs in the direction 
5G < 0, if the environmental change 5E is positive. Considering that evolution occurs 
through the projected direction in 5E, it is natural to assume 7 ^ = 7 ^ (although this 
might be a crude approximation). Under the linear conditions of interest, the change in p is 
proportional to 5E or 5G, with 6fi{6E,6G) = a{6E + 6G). Note that, again, the direction 
of SG is opposite to SE. Thus, we obtain, 

6X,{6E, 6G) = Lji{6fii6E, 6G) - ^i{6E + 6G)) = hp(hU, 6G) Y “ t/«)- (4) 

i i 

Then, over the course of evolution SG = 0 to 6G{m), under a given environmental condition 

E, 

_ 6Xj{6E,6G{m)) _ 6i^{6E,6G{m)) 

SXf'^^ ~ 5Xj{5E,Q) (5p((5U,0) ‘ 

In other words, all the expression changes are proportional, as conhrmed in the present 
simulations. To check the validity of the theory, we compared the proportion coefficient 
in the expression change (LHS of eq.5) with the change in growth rate (RHS) numerically 
through the course of the evolution simulation. As shown in Fig. 3, the relationship of eq. 
5 holds rather well. Note that if there is deviation from 7 ® = 'yf over i, the proportionality 
over all genes will deviate. In other words, the deviation from SXj{6E,6G) oc SXj{6E,0) 
across genes i reflects the deviation between 7 ^ and 7 ^. 

The relationship in the variances Vip{i) and Vg{i) is considered in a similar manner. 
Consider that the fluctuation in SE and SG induce fluctuation in each expression i, induced 
by either noise or genetic variation. This fluctuation induces variation in the growth rate, 
according to h/i = aSE or aSG, so that 

((iVUT)"> = LjXl- 7^/a))^ (6) 

i 

where hT is either SE or SG, i.e., phenotypic change induced by variation in the environment 
(i.e., noise) or by genetic change (e.g., by mutation), and (■ ■ ■) is the average over the 
distribution induced by the phenotypic noise or genetic variation. The variance Up(j) and 
Vg{j) are {{SXj{SE)Y) and (((5Xj((5G))^), respectively, so that 


u,(j) _ {Sii{SEf) _ u,(p) 


(7) 


Vg{j) {S^,{SGf) U,(/i) • 

Thus, the ratio of the two variances takes on the same value independent of j, which is 
determined by the ratio of variances in growth rate fluctuations induced by noise to those 
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induced by genetic variation. This relationship was again conhrmed in our simulated evolu¬ 
tion model (see Fig. S3) (see also 1^ for an alternative derivation of the common propor¬ 
tionality between Vip{i)/Vg{i) assuming the common error of catastrophe in the phenotype 
distribution). 

The above theoretical interpretations on both the proportionality in responses and in 
variances suggest the importance of changes in the growth rate. Since the growth rate glob¬ 
ally governs all of the concentrations through dilution, its dominance over each component 
is a reasonable assumption. To further evaluate the relationship between environmental 
and evolutionary dynamics, however, we need to also assume that evolution progresses 
as to assimilate environmental change, as Waddington proposed [l^. In our theory, this 
genetic assimilation is formulated by the introduction of the variable G that has a sim¬ 
ilar effect with the environment (or compensates for the environmental stress), so that 
dFi{xj)/dE ^ dFi{xj)/dG. 


Experimental Verification 

The theoretical argument and the simulation results demonstrated the existence of a 
common proportion coefficient r(m) for across components i, and its proportionality 

to the growth rate recovery To verify this relationship, we analyzed the 


time-series transcriptome data o 


conditions of ethanol stress 


19 


Dtained in an experimental evolution study of E. coli under 


20l |. In this experiment, after cultivation of approximately 


1,000 generations (2,500 hours) under 5% ethanol stress, 6 independent ethanol-tolerant 
strains were obtained, which exhibited an approximately 2-fold increase in specihc growth 
rates in comparison to the ancestor. For all independent culture series, mRNA samples were 
extracted from approximately 10® cells at 6 different time points, and the absolute expression 
levels were quantihed by using microarray analysis. All mRNA samples were obtained from 
the cells in exponential growth phase, which means that the changes in cellular state over 
the time scale of several generations were negligible, and each expression level represented 
cells in a steady-growth state (see |20| for details of materials and methods). 

Using the time-series expression data of bacterial adaptive evolution, we analyzed the 
common proportionality in expression changes. The environmental response of the Ath gene 
is dehned by the log-transformed ratio of the expression level of the Ath gene ob¬ 
tained 24 hours after exposure to the stress condition to that obtained under the no-stress 
condition. Similarly, the evolutionary response at n hours after the exposure to the stress 
5Xf-( n) is dehned by the log-transformed ratio of the expression level at n hours to that of 
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the non-stress condition. We found a common trend between the environmental and genetic 
responses over all genes, as shown in Fig. 6(a). Furthermore, we also found that the pro¬ 
portion coefficient r{n) for is roughly proportional to the growth recovery 

ratio as shown in Fig. 6(b), where ^/i‘^®"'(n) and are the growth rate 

differences of n hours and 24 hours after the exposure to stress, respectively. The results 
demonstrated that the evolutionary dynamics with growth recovery were accompanied by 
gene expression changes to eliminate the phenotypic changes introduced by the new envi¬ 
ronment, and agreed well with the simulation results of the simple cell model shown in Fig. 
2 as well as the theoretical argument presented above. 

Furthermore, there is some indirect experimental support for the proposed relationships 
between the variances. Stearns and colleagues measured the isogenic variance Vip of hve life- 
history traits (such as body weight, lifespan, etc.) in Drosophila melanogaster, as well as 
the genetic variance Vg between different genetic lines observed during laboratory evolution 
to increase the traits, and observed proportionality between the two 2l|]. The correlation 
between the isogenic variances in trait expression and variance due to mutation (but without 
selection) was measured across a few thousandgenes in Saccharomyces cerevisiae. Correla¬ 


tion between the two variances was observed 


22 l |. whereas proportionality was not so clear. 


This is possibly because evolution without selection was applied in the experiment, and 
therefore only the variance resulting from random mutation was measured. 


III. DISCUSSION 

We have shown proportionality in the change in the concentrations of most intra-cellular 
components as a result of adaptive evolution, which was conhrmed in simulated evolution 
of catalytic reaction network models of cells, laboratory experiments of bacterial evolution, 
and phenomenological theory. As the theoretical argument, albeit phenomenological, is 
rather general, we expect that the observed relationships obtained from the simulation and 
laboratory experiments represent a general phenomenon, independent of the specihc models 
or organisms considered. This proportionality across thousands of components implies that 
there is a strong constraint in phenotypic evolution. In particular, the expression of different 
components cannot evolve independently, but rather change together, for the most part, 
along a one-dimensional path provided by eq.(5). Phenotypic change in adaptive evolution 
under a fixed environmental condition is highly constrained; thus, we here quantitatively 
formulate the general restriction or feasibility of the direction of phenotype changes in future 
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evolution. 

Our result also implies that the changes in the concentrations of most components that 
are induced by the environmental change become relaxed through the evolutionary process. 
This suggests a phenomenon of strong homeostasis, that is, a tendency to restore to the 
original, adapted, intra-cellular states, via genetic change. In some sense, this homoeostasis 
is similar to the Le Chatelier principle in thermodynamics, in that changes introduced by 
external perturbations are relaxed by subsequent temporal evolution. 

There could be a huge variety of genetic changes that yield the phenotypic changes re¬ 
quired for adaptation. Indeed, in our simulations, there were a variety of network structures 
that could achieve phenotypic adaptation. When the simulation was run again with a differ¬ 
ent seed of random numbers for mutations, the resulting network (i.e., genotypes) was differ¬ 
ent in each run, but the change in concentrations (phenotypes) followed the proportionality 
given by eq.(5), independently of the specihc genetic changes occurring during evolution. 
Furthermore, in bacterial evolution experiments, the results from different strains tended 
to follow the same proportionality law described by eq.(5). It is interesting to note that 
such correlated change in expression levels by genetic changes is also suggested in several 
experiments 6|, l7|]. It will be important to further conhrm the relationship eq.(5) in more 
laboratory evolution experiments, and to also unveil the underlying genotype-phenotype 
map that achieves the common, restricted change in expression levels observed in the exper¬ 
imental data. 

We have also found proportionality in the fluctuations in expression levels across compo¬ 
nents. As expected from Fisher’s fundamental theorem of natural selection 23|, the higher 
the genetic variance, the higher the evolutionary rate. Hence, the proportionality between 
Vipii) oc Vg{i) suggests that a higher isogenic variance of a given expression level due to 
noise would be accompanied by a higher rate in the change in the expression level due to 
evolution. Hence, our results suggest that the direction of evolutionary change in phenotypic 
space is likely to be predetermined by the isogenic variance of expression level due to noise. 

According to our theoretical framework, the responses and fluctuations in expression 
levels are represented by the macroscopic growth rate and its fluctuation. Therefore, the 
relationship between the response and fluctuations, analogous to thermodynamics, is rep¬ 
resented by the landscape of the growth rate as a function of phenotype (expression level) 
and the environment, in contrast to the established htness landscape represented in genetic 
space proposed by Sewall Wright j^. We hope that the present study will provide a basis 
for the development of a future macroscopic theory for phenotypic evolution. 
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IV. METHODS: MODEL SIMULATIONS 


The cellular state can be represented by a set of numbers (iVi, iV 2 , • • • , Nk), where iVj is 
the number of molecules of the chemical species i with i ranging from i = 1 io K. For the 
internal chemical reaction dynamics, we chose a catalytic network among these k chemical 
species, where each reaction from some chemical i to some other chemical j is assumed to 
be catalyzed by a third chemical i, i.e., {i + £ ^ j + £). A catalytic network is chosen 
randomly such that the probability that any two chemicals, i and j, are connected is given 
by the connection rate p. Some resources (nutrients) are supplied from the environment 
by transportation through the cell membrane with the aid of some other chemicals that 
are named ‘transporters’. The concentrations of nutrient chemicals in the environment are 
kept constant, and they have no catalytic activity in order to prevent the occurrence of 
catalytic reactions in the environment. Through the catalytic reactions, these nutrients are 
transformed into other chemicals, including the transporters. Here, we assume that all of 
the K chemical species are necessary for cell division. Thus, cell division is assumed to occur 
when the minimum number of species exceeds a threshold M, i.e., min^^ iVj > M (in all 
analyses M is set to unity). Chosen randomly, the parent cell’s molecules are evenly split 
among the two daughter cells. In our numerical simulations, we randomly picked up a pair 
of molecules in a cell and transformed them according to the reaction network. In the same 
way, transportation through the membrane was also computed by randomly choosing from 
molecules within the cell and from nutrients in the environment. The parameters were set 
as A = 1000 and p = 0.01. 

We studied the evolution of the replication dynamics by generating slightly modihed 
networks and selecting those that grew faster. First, n parent cells were generated, and 
the connecting paths of catalytic networks were chosen randomly with connection rate p. 
From each of the n parent cells, L mutant cells were generated by randomly replacing mpK"^ 
reaction paths, where pK'^ is the total number of reactions and m is the mutation rate per 
reaction per generation. Then, reaction dynamics were simulated for each of the nL cells 
to determine the rate of growth of each cell; that is, the inverse of the time required for 
division. Within the cell population, n cells with faster growth rates were selected to be the 
parent cells of the next generation, from which nL mutant cells were again generated in the 
same manner. Throughout the simulation, the parameters were set as n = 1000 and L = 5. 
In the simulations shown in Fig. 1, the mutation rate m was set to 1 x 10“^. 

The environmental change is given by changing the nutrient concentration ratio in the en- 
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vironment. In the evolutionary simulation shown in Fig. 1, there are two nutrient chemicals, 
each associated with one transporter chemical. Before adding the new environmental con¬ 
dition (generation<0), the concentrations of these two nutrients in the environment (ci,C 2 ) 
were set to (0.5, 0.5), while after the environmental change (generation>0), they were set 
to (0.9, 0.1). In the result shown in Fig. SI, to add a variety of environmental changes, we 
randomly selected a nutrient chemical and a transporter chemical for this nutrient among K 
total chemical species. Then, the concentrations of the new nutrient cnew and the original 
nutrients were set to (ci, C 2 , CTVeto) = (0.45,0.45,0.1). We iterated the random addition of a 
nutrient 10^ times to obtain the average environmental response 
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Figure Captions 


Figure 1 

Growth rate (fitness) increase over generations from generation —3000 to 0. At generation 
0 , the nutrient concentrations in the environment were changed. As a result, the growth 
rate drastically decreased at generation 0, and later recovered by the evolutionary dynamics 
(see Methods for details). 

Figure 2 

Response by environmental change versus response by evolution. Relationship between the 
environmental response and genetic response (a), (b), and (c) show the 

plots for m = 5, 10, and 50, respectively. The black solid lines are y = x for reference. 

Figure 3 

The relationship between growth recovery rate and the proportion coeffi¬ 

cient r{m). The proportion coefficient r(m) was obtained by using the least-squares method 
for the relationship of SXf"''" and (5Xf®"'(m) for m = 1 200. The black solid line is y = x 

for reference. 

Figure 4 

The relationship between Vip{i) and Vg{i). The variances were computed by using the 
network and environment at generation 0 (before the environmental change) shown in Fig. 
SI with various mutation rates. Vip{i) and Vg{i) were calculated based on the simulation 
results of randomly generated 10® networks. The solid line is y = x for reference. 
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Figure 5 

Proportionality relationships among genetic/non-genetic fluctuations and responses hold 
over all components. The arrows indicate the proportional relationships. 


Figure 6 

Response by environmental change versus response by evolution in E. coli adaption to 
ethanol stress, (a) Relationship between environmental response and genetic re¬ 
sponse n) for n = 2496 as a representative example. and were 

calculated by the log-transformed expression ratio between before and 24 hours after and 
n hours after exposure to ethanol stress. The blue line is obtained by least-squares fitting, 
while the black line is y = x for reference, (b) Relationship between growth recovery rate 
^^Gen^^Env proportion coefficient r(n). The proportion coefficient r{n) was 

obtained by using the least-squares method for the relationship of 5X^^'" and 5Xf'^^{n) 
for n = 384, 744, 1224, 1824, and 2496 hours. The growth recovery rate 
was calculated based on the experimental measurements (see 20(] for details). Among the 6 
independent culture lines in 20(], the results of 5 culture lines without genome duplication 
are plotted. The black line is y = x for reference. 
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(a) 




Figure 7: Supplementary Figure 1. The relationships of to (a) Vip and (b) Vg. 
The variances are computed by using the network and environment at generation 
0 (before the environmental change). For the details of environmental changes to 
calculate see Simulation methods. 
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Figure 8: Supplementary Figure 2. The relationship between Vip and Vg after the 
environmental change. Vip and Vg before the environmental change (generation 0), 
immediately after the environmental change (generation 1), and after the adaptive 
evolution (generation 20) are plotted. After the environmental change, both Vip and 
Vg increase , and then recover the original levels after 20 generations of the evolution. 
The solid line \s y = x for reference. 
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Figure 9: Supplementary Figure 3. The relationship between the ratio of Vg{j)/V^p{j) 
and Vg{^) /Vip{ij). The variance ratio Vg{j)/Vip{j) is calculated by the least square 
method for all components. The data points are obtained with m = 10“° x 2^ for 
i = , 8 . 
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